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Background: Neutron star and supernova matter at densities just below the nuclear matter 
saturation density is expected to form a lattice of exotic shapes. These so-called nuclear pasta 
phases are caused by Coulomb frustration. Their elastic and transport properties are believed to 
play an important role for thermal and magnetic field evolution, rotation and oscillation of neutron 
stars. Furthermore, they can impact neutrino opacities in core-collapse supernovae. Purpose: In 
this work, we present proof-of-principle 3D Skyrme Hartree-Fock (SHF) simulations of nuclear pasta 
with the Multi-resolution ADaptive Numerical Environment for Scientific Simulations (MADNESS). 

Methods: We perform benchmark studies of ^°®Pb and nuclear ground states and calcu¬ 
late binding energies via 3D SHE simulations. Results are compared with experimentally measured 
binding energies as well as with theoretically predicted values from an established SHE code. The 
nuclear pasta simulation is initialized in the so-called waffle geometry as obtained by the Indiana 
University Molecular Dynamics (lUMD) code. The size of the unit cell is 24 fm with an average den¬ 
sity of about p — 0.05 fm“^, proton fraction of = 0.3 and temperature of T = 0 MeV. Results: 

Our calculations reproduce the binding energies and shapes of light and heavy nuclei with different 
geometries. Eor the pasta simulation, we find that the final geometry is very similar to the initial 
waffle state. We compare calculations with and without spin-orbit forces. We find that while subtle 
differences are present, the pasta phase remains in the waffle geometry. Conclusions: Within the 
MADNESS framework, we can successfully perform calculations of inhomogeneous nuclear matter. 

By using pasta configurations from lUMD it is possible to explore different geometries and test the 
impact of self-consistent calculations on the latter. 


I. INTRODUCTION 


In the high-density environment of neutron star 
interiors and core-collapse supernovae, nuclear matter 
is expected to assume a variety of exotic shapes at 
the liquid-gas phase transition. The different con¬ 
figurations are created by an interplay between the 
repulsive long-range Coulomb force and a short-range 
attractive nuclear force m\- At densities of around 
p ^ O.OI fm“^, spherical clusters of nuclear matter form 
a lattice surrounded by electron gas and neutron liquid. 
In a simple picture, with increasing density, the spheres 
merge into tubes that eventually transform into plates. 
As the density increases further, nuclear matter and 
neutron matter switch their roles resulting in tubes of 
neutron liquid and, at higher densities, bubbles enclosed 
by nuclear matter. For p ^ 0.12 fm“^, the neutron 
star interior is composed of a homogeneous mixture 
of neutrons, protons and electrons. In addition to the 
above shape sequence, many non-trivial geometries can 
be present and their similarity with different types of 
pasta (e.g. spheres = gnocci, tubes = spaghetti, planes 
= lasagna) lead to the terminology nuclear pasta phases. 
The relevant region for nuclear pasta in neutron stars lies 
between the outer core and the inner crust. Although 
the radial width of this region is only several hundred 
meters (in comparison to a neutron star radius of about 
10 km) its thermal and deformational properties can 
impact neutron star cooling mi], oscillations [6], spin 
[7] and magnetic field evolution |8l |9]. Understanding 


the physical characteristics of nuclear pasta is therefore 
an important step towards a correct interpretation of 
neutron star observables in connection with nuclear mat¬ 
ter properties and equation of state. For core-collapse 
supernovae (CCSN), pasta phases can form in the 
collapsing stellar iron core and during the post-bounce 
phase in the proto-neutron star m- The latter is the 
hot and compressed stellar core which is formed during 
the CCSN and left behind after the explosion. Neutrinos 
that diffuse from the proto-neutron star interior play a 
crucial role for the CCSN explosion mechanism mm- 
The knowledge of the neutrino mean free path in hot 
nuclear matter — that can be modified by nuclear pasta 
[iSfiTl] — is very important in numerical CCSN studies. 
In addition, pasta phases could have an impact on nucle¬ 
osynthesis in CCSN and neutron star binary mergers m- 

Different approaches are taken to study pasta phases. 
These include calculations in the liquid-drop model 
in na nil, Thomas-Fermi and Wigner-Seitz cell ap¬ 
proximations muHH], molecular dynamics (MD) and 
quantum molecular dynamics [22H24] studies, static 
Hartree-Fock [25H28] and time-dependent Hartree-Fock 
simulations m EO]. Studies are usually performed 
in the so-called unit cell filled with neutrons, protons 
and electrons with specific symmetry assumptions and 
boundary conditions. The pasta matter is then described 
as an infinite lattice of unit cells. By studying different 
configurations in the latter and comparing their total 
energies the ground state can be identified as the config¬ 
uration with the lowest energy. For numerical studies, it 
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is important to note the non-trivial role of the simulation 
volume. As only periodic geometries that fit into the 
unit cell can be explored, the size of the simulation cell 
must be sufficiently large to at least contain one period 
of the pasta structure. Even if the latter is fulfilled, 
effects of the finite volume such as dependence on the 
simulation space geometry [31] and numerical shell 
effects m can appear. As a consequence, the simulation 
volume has to be maximized to ensure that finite size 
effects are minimal. However, this usually comes with a 
significant increase in computational costs. 

The advantage of MD studies lies in their ability 
to simulate large systems where the length of the 
simulation space is several hundred fm [22] 132] [33] and 
therefore exceeds the size of a unit cell. This allows 
to study pasta structures that are less bound to the 
geometry and boundary conditions of the volume. 
Furthermore, it is possible to explore bulk properties 
such as electrical and thermal conductivities, shear 
and bulk moduli. However, although MD approaches 
can include quantum effects, the nucleon interaction is 
typically given by a schematic two-body potential. For 
self-consistent quantum calculations that account for 
Pauli blocking, spin-orbit forces and nucleon pairing, 
Skyrme Hartree-Fock-Bogolyubov (SHF(B)) and density 
functional theory (DFT) simulations are usually per¬ 
formed. A current drawback of these methods is their 
high computational cost. The consequence of the latter 
is that the system size that can be studied in SHF and 
DFT calculations is typically much smaller than for 
MD methods [28f[3Ql [34H36] . Therefore, for large-scale 
simulations, the applied numerical framework needs 
to be highly parallelized and should scale well. The 
Multi-resolution ADaptive Numerical Environment for 
Scientific Simulations (MADNESS) has been developed 
to efficiently solve these type of problems exactly and 
is designed to run on modern supercomputer facilities. 
With that, our aim is to apply MADNESS to perform 
large-scale 3D SHF simulation of nuclear pasta. In the 
current work we introduce our approach and perform 
benchmark studies of nuclear ground states. We then 
perform first pasta calculations using MADNESS and a 
converged MD simulation as starting point. 

The paper is structured as follows: We first give an 
overview of the Skyrme Hartree-Fock approach and the 
MADNESS computational environment in sections [TT| 
and |ml respectively. We continue with a description of 
how we solve the SHF equations in MADNESS in section 


in section El In section lyTRl we show our first nuclear 
pasta simulations with MADNESS and close with a 
summary in section ivnl 


IV and then present our results for nuclear ground states 


II. SKYRME HARTREE-FOCK 
CALCULATIONS 


In this section we provide a brief overview of the 
Skyrme Hartree-Fock method. For a more detailed 
discussion see e.g. [37H39] . Instead of solving the 
Schroedinger equation for the A-nucleon wavefunction 
and the corresponding Hamiltonian the Hartree-Fock 
approach approximates the ground state of the nuclear 
configuration by a single Slater determinant T. The lat¬ 
ter is formed by a complete orthonormal set of single¬ 
particle wavefunctions 'ipi{fi) whereas fi contains the spa¬ 
tial, spin, and iso-spin coordinates of the ith state. The 
Slater determinant is obtained via the variational princi¬ 
ple by the minimization of the energy expectation value: 

6E{^) =6{^\H\^) =0. (1) 

As a result, the many-body Schroedinger equation is 
turned into A single-body Schroedinger equations with 
the Hartree-Fock mean-field potential I/hf- For each 
single particle state the corresponding HF equation 
reads: 

H V’i(jO = Ei (2) 

where rrii is the nucleon mass and Ei is the energy of 
the single-particle state. For an iso-spin state q = n^p 
(n=neutrons, p=protons), the Hartree-Fock potential 
U}iF,q contains the following contributions: 

U}lF,q = Uq ,sky T Vgf,meff T + current + ^g,spin- 

(3) 

Furthermore, for protons, the Coulomb potential Uq and 
Coulomb exchange potential Uc,ex are added. The first 
term in eq. © is the Skyrme potential and can be ex¬ 
pressed as 


Ug,sky = bop - b'opq + biT - b^Tg - b2^p 
+ h'^^pg + bo'^^-^p“+^ - h'^‘^p“pg 

- /5““' {pI + pI) - &4V ■ J - ■ Jg. (4) 

Here, ct, bj and 6'- {j = 0...4) are constants specific to the 
Skyrme potential. They are fitted to reproduce known 
properties of finite nuclei and infinite nuclear matter. 
The nucleon number densities pq and kinetic densities 
Tq in eq.Q are given by: 

p? = FF( 5 ) 

is is 

Np = Nn = A — Z ^ 

where A is the mass number, Z the charge number of the 
nuclear configuration and s marks the spin of the state. 
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The divergence of the spin-orbit density J is determined 
by: 

^ E E X • (s' I I s)) • (6) 

i ss' 

where are the Pauli matrices. The Coulomb exchange 
potential Uc,e^ in eq.(|^ can be calculated via the so- 
called Slater approximation [40] : 

\ 1/3 

Ppi^) • (7) 

The Coulomb potential of the protons is given by: 

Ucp = f ds. (8) 

Furthermore, when studying nuclear pasta phases, we as¬ 
sume that the unit cell is charge neutral and contains neu¬ 
trons, protons and electrons. The latter are given via the 
so-called Jellium approximation and form a background 
of homogeneous negative density pj = —ZjV . As a con¬ 
sequence, in addition to proton Coulomb potential, we 
have to consider the interaction of the protons with the 
Coulomb potential of the Jellium: 

Ucj = [ ds. (9) 

J |r-Jl 

We can sum both contributions to a total Coulomb po¬ 
tential: 

Ucp + Ucj = J """ ^ ds, (10) 

where 

Uc = Ucp + Ucj, Pc(f) = Pp(f) + pj- (12) 

We will apply the Jellium approximation whenever 
studying volumes with periodic boundary conditions. 
The remaining components of the nucleon potential in 
eq.Q are a contribution that accounts for the effective 
nucleon mass: 

= -V • {bip - b[pq) V, (13) 


play m]. 

Since the single particle states in eq. <§ depend on 
the potentials and densities that in turn are derived from 
the wavefunctions, HF problems have to be solved itera¬ 
tively with an assumption about the initial single-partilce 
states e.g., harmonic oscillator states, 3D gaussians or 
plane waves. From these, we derive densities and poten¬ 
tials that are used in the SHF equations to determine a 
new set of updated states. The calculations are repeated 
until the solution for the wavefunctions is self-consistent. 
Due to the large number of states in nuclear pasta simu¬ 
lations, we require a numerical framework that is compu¬ 
tationally efficient and parallelized. We therefore apply 
the Multi-resolution ADaptive Numerical Environment 
for Scientific Simulations (MADNESS) which will be de¬ 
scribed in the next section. 


III. MULTI-RESOLUTION ADAPTIVE 
CALCULATIONS 

MADNESS is a numerical framework designed to effi¬ 
ciently solve problems involving integral and partial dif¬ 
ferential equations in many dimensions. Examples in¬ 
clude Hartree-Eock and density functional theory calcu¬ 
lations of chemistry and nuclear physics problems [E]- 
|48] with a recent application in studying finite nuclei via 
solving the HEB equations [42]. Operations in MAD¬ 
NESS are highly parallelized via a combination of MPI 
and ptbreads parallel computing. 

In MADNESS, functions and operators are described by 
adaptive pseudo-spectral approximations that are based 
on a multi-wavelet basis. The latter is given by discon¬ 
tinuous Alpert’s multi-wavelets ns EOl with Legendre 
polynomials being applied as scaling functions. Both, 
scaling functions and multi-wavelets, have disjoint sup¬ 
port and are efficient in describing discontinuities and 
regions with high curvature. Eurt her more, with each op¬ 
erator and function having its own adaptive structure of 
refinement, the user can achieve a defined finite but guar¬ 
anteed precision. In the following, we will briefly describe 
the multi-resolution approach in MADNESS whereas de¬ 
tails can be found in e.g. lEiisoiini]. 

MADNESS projects functions and operators from the 
user space with a defined width onto a solution interval 
[0,1]. Here, k orthonormal Legendre scaling functions 


^C,ex(^) — 6 I 

\ TT 


and the spin-orbit potential: 

Uq^so = Wq • (a X V) , W, = ( 64 Vp + b'^Vpq) (14) 

In this work we are focusing on time-independent HE cal¬ 
culations of even-A and even-even nuclei as well as pasta 
phases. Therefore, we do not include the current and 
spin operators I/g,current and I/g,spin, respectively. These 
contribute to the Hamiltonian only in case of odd-A and 
odd-odd nuclei, and when dynamical effects come into 


V 2 i + lPi{2x — 1 ) for 0 < X < 1 
0 otherwise ^ ' 

i = 0 ,..., k — 1 

can be defined. They are the /th Legendre polynomi¬ 
als Pi{x) shifted to [0,1] and normalized. The solution 
interval is repeatedly cut in half. At level n, there are 
2^ boxes of size 2^~^. The functions 0i(x) are scaled 
to level n and translated to each subinterval I with size 
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[2 2 + 1)] where they are given by: 


4>l{x) = V¥^ct>i{2'^x-l), (16) 

i = 0, — 1, / = 0, — 1 

The scaling functions are orthonormal on the interval 
2~^{l +1)] and span the sub-spaces which form 
a ladder: 


c V*' c V 2 ^ ... c c ... (17) 


Due to this relation, scaling functions at level n can be 
derived by scaling functions at level n + 1 by the two-scale 
relationship [52] . In the reconstructed form, a function / 
that is smooth at level n, can be represented by scaling 
functions and coefficients s'^i as: 


2^-1 k-1 

r{x) = 

Z=0 i=0 
p 2 ^(Z+1) 


^il = 


/ 

J2-^l 


4>u{x) f{x)dx. 


(18) 

(19) 


ID to 3D, functions are given by tensor products of 
multi-wavelets and scaling functions are given in the non¬ 
standard form. Adaptive refinement is performed locally 
if the local error is above a truncation threshold e. In 
ID, it is accomplished by truncation of small wavelet co¬ 
efficients whereas MADNESS offers different truncation 
criteria. For an accurate representation of functions and 
their derivatives, which we are interested in, coefficients 
for level n and sub-interval I are neglected when: 

IMrih = E ^ ^ min(l,2-V), (25) 

i 

where L is the minimum width of the simulation volume 
and e is the desired precision. In MADNESS, Green’s 
functions are represented via low-separation rank expan¬ 
sion in terms of Gaussians. For example the Yukawa 
kernel is: 


^ — k r 


M 


= E 




-Pi, 


0-^2,' 


o-PS,^ 


m=l 


+ o{- 

^r 


(()■ 


(26) 


The complementary subspace to in is with: 

= (20) 

It is spanned by multi-wavelets on the interval 

[2“^/, 2“^(/-|-l)] that are obtained by dilation and trans¬ 
lation of ^i'. 

i>l{x) = V^i>i{2^x-l), ( 21 ) 

i = 0,...,fc-l, ^ = 0,...,2”-l 

which, in turn can be derived from the multi-scaling func¬ 
tions by the two-scale relations. Alpert’s wavelets are or¬ 
thonormal within and between scales. Since can be 
decomposed into: 

y* = © W'l © ... © 1, (22) 

the function can be given as a sum over scaling func¬ 
tions at the coarsest level and wavelets at finer length- 
scales: 

k-1 / n-1 2^-1 \ 

/"(y = E ( 4>i{x) + E E ’ (^^) 

i=0 \ m=0 Z=0 / 

p2 ”^(Z+1) 

= / fix)ij^dx ( 24 ) 

J2-^l 

This is the so-called compressed form. The reconstructed 
representation and compressed representation are two 
equivalent forms of /^. For some numerical operations it 
is better to use the scaling function representation while 
for others (e.g. inner product of functions) the wavelet 
form is more efficient. The transformation between both 
representations of is an orthogonal transformation, 
it is therefore numerically stable and fast. Going from 


This reduces a 3D convolution to a set of uncoupled ID 
convolutions with M depending on the user-determined 
precision e. Transformation matrices with respect to the 
multi-wavelets are pre-computed which allows a fast com¬ 
putation of the convolution. 

IV. SOLVING THE SKYRME HARTREE-FOCK 
EQUATIONS 

Our general strategy is to rewrite a given differential 
problem into an integral form and solve it iteratively 
via convolutions with Green’s functions. Gorrespond- 
ingly, we rearrange each HF equation in eq. into their 
Lippmann-Schwinger form: 

{a~^A += UMF,q{r)'ipiifl, (27) 

where a~^ = /2m. This can now be expressed as a 

convolution with the Green’s function for the bound state 
Helmholtz (BSH) equation: 

= -OL GbSH,! ^ {UuF,q (28) 

/ oo 

GBSH,i(G 5 ) (GHF,q (5)'0i(s)) ds (29) 

-CX3 

GBSH,i(r,s) = , |f '"-*1, k = (30) 

47r|r — ^ 

Similarly, the Goulomb potential Uc for a given total 
charge density pc and 

MJc = -4:Trpc (31) 

is given by the convolution with the corresponding 
Green’s function Gc{r^s) = l/\f— 
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To smooth out possible numerical noise in first and sec¬ 
ond derivative terms, we apply Gaussian smoothing. Due 
to its high resolution and adaptive refinement, MAD¬ 
NESS resolves small discontinuities which could then 
propagate from e.g. the Skyrme potential into the up¬ 
dated wavefunctions. If not damped, the noise can am¬ 
plify with each iteration. As a consequence, we convo¬ 
lute kinetic densities Tq and density laplacians with 
Gaussians of the following form: 

/(r) = \ (32) 

where typically 0.25 fm is chosen for a. Smoothing is 
applied for the initial iterations and removed once the 
configuration starts to converge. 

Our code is based on a previous algorithm to solve 
the Lippman-Schwinger equations for HE problems with 
spin-orbit potential [44]. The iterations are therefore 
performed in a similar fashion. We start out with a 
set of single particle states at iteration n = 0 - 
These states are ortho-normalized using the a LAPACK 
hermitian eigensolver for the generalized eigen-system 
problem [44| 

HC = SCE, Hij = J (33) 

Sij = I (34) 

for C and E. The new states with energies Ei are 
obtained via 

= (35) 

We then apply the potential operator I/HF,q and perform 
the convolution with Gbsh,!- This gives us a new set of 
states at iteration n + 1: 

0n+l ^ GBSH,i ^ (GHF,q ) (36) 

= -a{-A-aEi)-^ C/nF.q </>”, 

UuF,q = Uq-V ■ {hip - h\pq) V + iWq • (a X V). (37) 

To determine the convergence of the states we calculate 
the maximal L2-norm of the wave-function difference be- 


tween two iterations: 


dij) = 

('5V’i)max = TXiax{6lpo,...,6'lpNg}, 

(38) 

5ipi = 


(39) 


If is smaller or equal to a given desired precision e, 
the iterations are considered as converged. Otherwise, 
the new single-particle states are calculated as averages 
of the old and new wavefunctions: 

(40) 


and the next iteration is performed. The averaging is a 
usual technique in HE calculations to stabilize the iter¬ 
ations. In our simulations we typically use y = 0.4. A 
larger value of y leads to a faster convergence but might 
also allow the development of instabilities. Note that 
there are different iteration routines as well as conver¬ 
gence criteria I271I53H55] which might be more suitable 
for SHE calculations and will be tested in the future. 
Eor the present study we apply the same iteration steps 
and check for convergence as has been done in previous 
MADNESS studies [44] . 


V. NUCLEAR GROUND STATES 

Before we apply our code to study nuclear pasta 
phases, we want to test its performance and accuracy 
by simulating the known ground states of several nuclei 

- ^^^Pb, and The first, is a doubly 

magic nucleus with a well-known binding energy and 
is therefore a good first benchmark test of our code. 
Similarly, ^^^Pb is also doubly magic but with 13 times 
more nucleons than Einally, is a deformed 

nucleus and thereby a good test case for our code to 
find the nuclear ground state throughout several shape 
changes. 

We calculate nuclei in (I) a large simulation space 
with free boundary conditions (be) and no Jellium, 
and (H) a small box with periodic be including the 
Jellium approximation. While the first is suitable 
for comparisons with experimental ground states, the 
second case has similar conditions as our nuclear pasta 
simulations. Eurthermore, in addition to experimental 
binding energies [56], we compare our results with the 
Skyrme HE code Sky3D [29] [^ [35l [37] . Eor the setup 
of simulation (I) in M-SHE, we choose a box width of 
L = 200 fm, while for (H) we apply L = 24 fm. The latter 
is used in Sky3D for simulations (I) and (H). Although 
L = 24 fm is a relatively low value, the nuclear densities 
are ^ 10“^ fm~^ at the edges of the simulation space 

- even for large nuclei. However, to cross-check that 
the small simulation volume does not alter the nuclear 
ground states for simulation (I), we perform calculations 
with L = 48 fm and compare the total energies to 
L = 24 fm. Einally, to test the dependence on resolution 
in Sky3D, we vary the cell size of the computational 
grid. Eor MADNESS, we decrease the final truncation 
threshold e and check the impact of Gaussian smoothing. 

The single-particle states with spin-up and spin-down 
are initialized as harmonic oscillator states: 

'Ipi,s{r^ = i = 0,...,Ns-l (41) 

where j, k and I are integers starting at 0, Ng is the 
number of states with spin-up or spin-down, and d = 
0.625 ^1/3. We use the Skyrme force SV-bas |57] and 
calculate the binding energy from the energy components 


+(1 - x)<^r 
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of the Skyrme density functional: 

Eo = l [bo (rl -b'oJ2pl (42) 

-El = y 6ip(f)r(f) - 6'i y] pq{r)Tq{f) dr, (43) 


Eo — — 


\ J - ^2 p{r)^p{r) dr, (44) 

Eo = ^Jhp^+\f^- b’, p^{E)Y1 pI ( 45 ) 


E4 = - 


j hi p{r) VJ{r) + h'^Y. Pi(^ (46) 


the kinetic and Coulomb energies: 

^kin = y] ^ y Tq{E) dr, ( 47 ) 

^C = IJ Uc{i^pp{i^ df. ( 48 ) 


For the latter, the total charge potential Uc and pc are 
given by eq.(12) for simulation (II) and Uc = Upp for (I). 
The total binding energy and binding energy per nucleon 
are then given by: 




total 


• ^kin + Ec + + ^2 + ^3 + ^4 (49) 


and F^bind = F^totai/A, respectively. As in Sky 3D, we do 
not consider energy contributions beyond the mean-field 
approximation, for example a center-of-mass correction 
m- The latter decreases with higher mass number as 
^ 1/A and should therefore have only a small contribu¬ 
tions in simulations of heavy nuclei and nuclear pasta. 
However, for light nuclei, the lack of a center-of-mass 
correction might lead to noticeable deviations from ex¬ 
perimentally obtained binding energies. Results of the 
simulations are given in tables 0-111 Here, for Sky 3D 
and the MADNESS SHF code (abbreviated as M-SHF), 
we give the box length L together with the simulation 
setup (I) or (H). The resolution of the simulation is given 
by the grid cell size Ax for Sky3D and truncation thresh¬ 
old e for MADNESS. The binding energy and total bind¬ 
ing energy, Fbind and FtotaU respectively, are followed by 
the different energy components as described in eq.(42) - 
eq.i 


_ |). All energies are given in units of MeV. For the 

nuclear simulations with M-SHF, we typically start out 
with e = 10“^ and Gaussian smoothing using a = 0.25. 
When 


d^p rsu A X e (50) 

we decrease the truncation threshold to e = 10“^. At this 
point, we continue with three versions of the simulation. 


The first version still contains the Gaussian smoothing 
while we remove it in the second one. Both simulations 
are evolved until eq.(50) is fulfilled for the new e. At this 
point, the calculations are stopped. We can then com¬ 
pare the impact of the smoothing on the energies and the 
nuclear configuration. In a third simulation, we continue 
the simulation without smoothing and decrease e by fac¬ 
tors of 10 according to eq.(50) while e > 10“^. The simu¬ 
lations are run on the high performance computer center 
BigRed H at Indiana University and the EOS cluster at 
the Oak Ridge Leadership Gomputing Facility where we 
use nodes with 16 cores. As our code is not yet opti¬ 
mized for speed we do not give specific runtime numbers 
at this point but provide typical order-of-magnitude it¬ 
eration counts and computational times for each nucleus 
and pasta simulations. Note that the time-independent 
calculations with Sky3D are not MPI parallelized and 
run at BigRed H on one node with 32 cores. As Sky3D 
represents wavefunctions and performs calculations on a 
fixed grid, simulations times for a given number states 
scale roughly by a factor of eight when L is increased 
by a factor of two or Ax decreased to half its size. Due 
to the adaptive refinement, MADNESS simulation times 
depend mostly on the truncation threshold e (again as¬ 
suming a fixed number of wavefunctions). 


A. nucleus 

Our results for the nucleus are given in table 

for the Sky3D and M-SHF simulations. We first discuss 
the results for setup (I). Figure [^(a) shows the num¬ 
ber density profiles of the converged nucleus for M- 
SHF with e = 10“^ along the three axes and —10 fm 
< x^y^z < 10 fm. Due to the spherical shape of the 
nucleus, the profiles overlap exactly. The kinetic density 
and laplacian of the total density p are shown in Fig.[^(b) 
along the x-axis, together with the corresponding density 
profile. Despite the absence of Gaussian smoothing, the 
profiles show no discontinuities or irregularities and, due 
to the spherical symmetry of the nucleus, are identical 
along all axes. In table [I[ we find a clear difference be¬ 
tween energies for simulations with Gaussian smoothing 
(marked with a in the resolution column) and without, 
whereas all components are affected. The difference in 
F^totai between a simulation using smoothing and without 
is |AF;totai| ^ 1-458 X 10“^ |F;totai| for e = 10“^, whereas 
this value will depend on the size of cr. On the other 
hand, when comparing simulations without smoothing 
but with different truncation thresholds, e = 10“^ and 
e = 10“^, we do not find any noticeable differences for 
up to 7 decimals in Ftotai (b decimals in the table). 

For Sky3D, changing the box size from F = 48 fm to 
F = 24 fm leads to a difference in total energy of only 
I^F^totail ^ 1-20 X 10“^ |Ftotai|- This is smaller than the 
change due to a decrease in cell size from Ax = 1 fm to 
0.5 fm which results in |AFtotai| ^ 3.86 x 10“^ |Ftotai|- A 
further reduction to Ax = 0.25 fm has a negligible effect 
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FIG. 1: (a) Number density profiles of nucleus ground 
state, as obtained with our MADNESS code, taken along the 
X, y, and z-axis. The profiles show that the resulting nucleus 
is spherically symmetric, (b) Number density, kinetic density 
T and laplacian of the number density Ap along the x-axis. 


which implies that, at least for Ax = 0.5 fm is suf¬ 
ficient to capture the correct energy values. 

Regarding their sensitivity to resolution, results for sim¬ 
ulations (II) are very similar to the ones of setup (I). 
From table [T| we see again a difference in the total en¬ 
ergy for Sky3D simulations when changing the cell size 
from Ax = 1 fm to 0.5 fm, and no effects for a further de¬ 
crease to Ax = 0.25 fm. Similarly, for M-SHF, Gaussian 
smoothing leads to |AF^totai| ^ 3.51 x 10 ^ |£'totai|, while 
a decrease in truncation threshold does not have any vis¬ 
ible effects. However, for all simulations, the absolute 
values of the total energy are higher than in setups with 
free be. This is due to the smaller Coulomb energy as a 
consequence of the Jellium. The latter reduces the total 
electric charge density and thereby the Coulomb poten¬ 
tial. Although M-SHF requires less than hundred itera¬ 
tions until convergence while Sky3D uses several hundred 
steps, both codes are very fast and require less than one 
hour on one node. Figureshows the maximum error 
and the binding energy per particle for M-SHF as they 
evolve with iterations for setup (I) and a final truncation 
threshold of e = 10“^. The vertical dashed lines mark 



Iterations 

FIG. 2: Evolution of the maximum error and binding energy 
with iterations of the calculation. 


the reduction of the truncation threshold to the new val¬ 
ues as given in the figure. The large jump in Sip and 
^bind at iteration ^ 30 is due to the removal of Gaussian 
smoothing. After that, we see that the binding energy 
does not change much until the calculation is converged 
and Sip becomes constant. The MADNESS calculations 
of and the Sky3D results are in good agreement with 
each other for simulations (I) and (H) and differ by only 
lAE’totail ^ 0.002 MeV when using the highest discussed 
resolutions. The large deviation from the experimental 
binding energy of F^exp ^ —7.976 MeV [56] originates in 
the applied Skyrme force Sv-bas and, as previously men¬ 
tioned, could be partially attributed to the absence of the 
center-of-mass correction. 


B. 2°®Pb nucleus 

Next, we discuss simulations of the ^^^Pb nucleus 
with Sky3D and M-SHF. The resulting energies are 
summarized in table an where the structure of the table 
is as for Interestingly, for Sky3D, the difference 

in the total energies of ^^^Pb between L = 24 fm and 
L = 48fm in simulation (I) is the same as for namely 
lAE’totail ^ 0.014 MeV. In comparison to the total en¬ 
ergy of ^^^Pb it is of course only ^ 8.53 x 10“^|F^totai|- 
When increasing the resolution by setting Ax = 0.5 fm, 
the total energy changes by about ^ 0.058 MeV or 
^ 3.56 X 10“^F^totai- As before, it seems that the change 
in energy due to the simulation volume is smaller than 
the one caused by an increase in resolution. Setting 
Ax = 0.25 fm results in |AF^totai| ^ 0.003 MeV which 
is negligible in comparison to the total energy. As for 
we can conclude that a resolution of Ax = 0.5 fm 
is sufficient to reproduce the ground state of ^^^Pb and 
a cell size of L = 24 fm does not lead to large finite-size 
effects. For ^^^Pb in simulation (H) and the 
calculations we will therefore only test L = 24 fm and 
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L [fm] 

sim. 

resol. 

Ebind 

Etotal 

Ekin 

Eo 

El 

E 2 

E 3 

E 4 

Ec 

Sky3D 

48 

(I) 

1.0 fm 

-7.290 

-116.643 

234.538 

-976.171 

12.689 

43.471 

556.049 

-0.747 

13.542 

Sky3D 

24 

(I) 

1.0 fm 

-7.291 

-116.657 

234.537 

-976.163 

12.689 

43.470 

556.044 

-0.747 

13.542 

Sky3D 

24 

(I) 

0.5 fm 

-7.288 

-116.612 

234.443 

-976.109 

12.687 

43.499 

556.086 

-0.752 

13.535 

Sky3D 

24 

(I) 

0.25 fm 

-7.288 

-116.613 

234.443 

-976.113 

12.687 

43.499 

556.088 

-0.752 

13.534 

M-SHF 

200 

(I)* 

10“'’ 

-7.396 

-118.336 

237.165 

-992.391 

12.801 

43.528 

567.719 

-0.765 

13.606 

M-SHF 

200 

(I) 

10“® 

-7.288 

-116.611 

234.444 

-976.114 

12.688 

43.499 

558.088 

-0.752 

13.535 

M-SHF 

200 

(I) 

10“'^ 

-7.288 

-116.611 

234.444 

-976.114 

12.688 

43.499 

556.088 

-0.752 

13.535 

Sky3D 

24 

(ii) 

1.0 fm 

-7.626 

- 122.010 

234.606 

-976.537 

12.697 

43.497 

556.312 

-0.747 

8.192 

Sky3D 

24 

(ii) 

0.5 fm 

-7.622 

-121.958 

234.504 

-976.439 

12.694 

43.522 

556.322 

-0.752 

8.191 

Sky3D 

24 

(II) 

0.25 fm 

-7.622 

-121.958 

234.504 

-976.438 

12.694 

43.522 

556.321 

-0.752 

8.191 

M-SHF 

24 

(II)* 

“DW 

-7.730 

-123.683 

237.222 

-992.699 

12.807 

43.549 

567.940 

-0.764 

8.262 

M-SHF 

24 

(II) 

10“® 

-7.622 

-121.957 

234.505 

-976.440 

12.694 

43.522 

556.321 

-0.752 

8.192 

M-SHF 

24 

(II) 

10“'^ 

-7.622 

-121.957 

234.505 

-976.440 

12.694 

43.522 

558.321 

-0.752 

8.192 


TABLE I: Parameters and energies for SkySD and M-SHF simulations of the nucleus. Simulations with free boundary 
conditions are marked by (I) while periodic boundary conditions with the jellium approximation are given by (II). The simulation 
box size is given by its length L. The resolution is defined as the grid cell size for SkySD and truncation threshold for M-SHF. 
Simulations that apply Gaussian smoothing are marked by a The bind ing energy per baryon Fbind, total energy Ftotai and 
different energy components: Fkin, Eq - E 4 , and Ec (see eq.(42) - eq.(48)) are given in MeV. Binding energies for simulations 
that were performed with the highest precisions are marked by bold font. 


Ax > 0.5 fm. As we will see, the energetic differences 
between simulations with Ax = 0.5 fm and Ax = 1.0 fm 
are very small. 

For M-SHF, we apply again Gaussian smoothing with 
a = 0.25 in the beginning of the simulation when we 
initialize the wavefunctions as harmonic oscillator states 
and truncate with e = 10“^. Once eq.(50) is fulfilled, 
we continue with three simulations as described in the 
previous section. Table [ll| shows that Gaussian smooth¬ 
ing affects again all energy terms, leading to a difference 
in F;totai of ^ 1.674 MeV or ^ 1.02 x lO^^lF’totail- In 
contrast, reducing e from 10“^ to 10“^ results in a 
small change of the total energy of only ^ 0.002 MeV or 
1.22 X 10 ^l^totail- This suggests that a final value of 
e = 10“^ is sufficient to reproduce the nuclear ground 
state energies. 


Figure |^a) shows the x, y and 2 : profiles of the 
initial and final total density p in simulation (I) for 
e = 10“^ and —20 fm < x^y^z < 20 fm. We can see 
that while the initial density distribution is slightly 
flatter in the 2 :-direction, the final shape is spherically 
symmetric. The x-profile of the total density is again 
shown in FigJ^b) together with the profiles for r and 
Ap. We also plot the corresponding final quantities for 
simulation (I) with e = 10“^ and Gaussian smoothing. 
Differences in p and r are not visible while the Ap 
profiles show small deviations around x ^ 0 fm. For 
better comparison, we zoom into the nucleus and 
compare the profiles again in Figj^ Small deviations in 
the oscillation amplitudes of all three quantities can be 
seen. However, despite these differences, the oscillation 
pattern are very similar. As a consequence, while the 
application of Gaussian smoothing impacts the energies, 
the effects on the shape of a nuclear configuration might 
be small. 


As for the sensitivity of the periodic simu¬ 

lations for ^^^Pb is similar to the free be stud¬ 
ies. For Sky3D, increasing the resolution by 

changing Ax from 1.0 fm to 0.5 fm, results in 
lAEtotail - 0.087 MeV - 4.03 x 10“^ |F;totai|. With 
lA^totail - 1-694 MeV - 7.84 x Gaussian 

smoothing in M-SHF has again a larger impact on the 
energies than reducing the truncation threshold which 
results in |AEtotai| ^ 0.002 MeV. 

In general, both codes agree well. The deviations 
between numerical results for F^bind in simulation (I) 
and the experimental value of ^ —7.867 MeV are due to 
the applied Skyrme force SV-bas. The higher energies 
in simulations with periodic boundary conditions are 
again due to the inclusion of Jellium. Although, the 
difference in total energy between both codes results in 
only ^ 0.8 MeV or ^ 0.004 MeV per baryon, it is notable 
and should be understood. As can be seen from tables 
|T] and |TI| the difference AF^totai between Sky3D and 
M-SHF seems to scale with the number of states and is 
most pronounced in the Eq and E^ terms which are both 
functions of the baryon densities and have large absolute 
values. Our tests of L, Ax and e did not reveal any 
sensitivities of the results that would be large enough to 
account for the energy difference. Furthermore, since we 
see the same behavior for simulations (I) and (II), we 
can assume that it is not an effect of specific boundary 
conditions. With that, more cross-checks have to be 
performed between both codes in the future. 

Due to the larger number of wavefunctions, calcu¬ 
lations of ^^^Pb take longer than for Simulations 

with Sky3D, required ca. 3.5 hours and about 1000 
iterations on one node for Ax = 1.0 fm. With M-SHF, 
simulations of ^^^Pb were performed on 8 CPU nodes. 
Gaussian smoothing has only a minor impact on the 
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L [fm] 

sim. 

resol. 

Ebind 

Etotal 

Ekin 

Eo 

El 

E 2 

E 3 E 4 Ec 

Sky3D 

48 

(I) 

1.0 fm 

-7.842 

-1631.027 

3920.428 

-17379.720 

285.346 

240.119 

10591.390 -86.064 798.492 

Sky3D 

24 

(I) 

1.0 fm 

-7.841 

-1631.013 

3920.767 

-17380.820 

285.356 

240.123 

10591.170 -86.059 798.489 

Sky3D 

24 

(I) 

0.5 fm 

-7.841 

-1630.955 

3920.585 

-17380.620 

285.356 

240.101 

10591.220 -86.067 798.465 

Sky3D 

24 

(I) 

0.25 fm 

-7.841 

-1630.958 

3920.587 

-17380.640 

285.357 

240.102 

10591.240 -86.067 798.462 

M-SHF 

200 

(I)* 

mw 

-7.845 

-1631.834 

3923.895 

-17397.536 

285.436 

239.842 

10603.944 -86.106 798.692 

M-SHF 

200 

(I) 

10 “® 

-7.837 

-1630.160 

3920.557 

-17376.751 

285.309 

240.030 

10588.270 -86.045 798.469 

M-SHF 

200 

(I) 

10 “'^ 

-7.837 

-1630.162 

3920.541 

-17376.675 

285.309 

240.021 

10588.219 -86.047 798.469 

Sky3D 

24 

(ii) 

1.0 fm 

-10.385 

-2160.078 

3940.437 

-17519.980 

289.937 

242.621 

10701.650 -87.118 272.412 

Sky3D 

24 

(ii) 

0.5 fm 

-10.385 

-2159.991 

3940.278 

-17520.010 

289.943 

242.615 

10701.890 -87.127 272.419 

M-SHF 

24 

(II)* 


-10.389 

-2160.857 

3943.894 

-17537.870 

290.025 

242.342 

10715.300 -87.162 272.616 

M-SHF 

24 

(II) 

10 “® 

-10.381 

-2159.163 

3940.599 

-17517.276 

289.908 

242.549 

10699.733 -87.100 272.424 

M-SHF 

24 

(II) 

lO”’^ 

-10.381 

-2159.161 

3940.577 

-17517.170 

289.907 

242.540 

10699.663 -87.101 272.424 


TABLE II: Parameters and energies for SkySD and M-SHF simulations of the nucleus. The table setup is the same as 

in tableThe experimental binding energy for ^°®Pb nucleus is Ebind ~ —7.867 MeV [56] . 



distance [fm] distance [fm] 


EIG. 3: (a) Number density profiles of the ^°®Pb along the x, y, and z-axis in the initial state (subscript i) and at convergence 
(subscript /). Subfigure (b) shows the number density, kinetic density r and laplacian of the number density Ap at iteration 


280. 


simulation time while setting the truncation threshold 
to lower values increases the latter. However, in all 
cases, calculations reach the ground state within several 
hours and < 300 iterations. The binding energy and 
log((5'0) as functions of iteration number for simulation 
(I) and final truncation threshold e = 10“^ are plotted 
in Fig.[^ The binding energy is calculated every 10 
iterations and the dashed vertical lines mark again the 
reduction of e. We can see three peaks in d'ljj which 
are caused by the decrease of the truncation threshold 
according to eq.(50). The first decrease is accompanied 
by the removal of Gaussian smoothing which also leads 
to a jump in the binding energy. However, except for 
the three discontinuities, the value of gradually 
decreases, which indicates that the initialization via 
harmonic oscillator states is a good initial guess and 
does not require any major shape changes of the nucleus. 


C. 238^ nucleus 

Finally, we discuss the nucleus. The density 

profiles for the initial and final states in simulation (I) 
are shown in FigJ^a) for —20 fm < x, p, ^ < 20 fm. As 
before, we start with harmonic oscillator states. The 
initial density distribution of is again squeezed in 
the z-direction while the final nucleus is elongated along 
the y-axis. The many required shape changes result in 
a longer convergence time as will be discussed at the 
end of this section. Results for different simulations 
with Sky3D and M-SHF are given in table |HI[ As 
previously mentioned, for Sky3d, we perform simulation 
with Ax = 1 fm and L = 24 fm. A smaller value of Ax 
and larger simulation space results in small changes in 
the energy contributions. As can be found in table |HI[ 
final truncation thresholds of e = 10“^ and e = 10“^ in 
M-SHF do not lead to large differences in the energies. 
Gaussian smoothing, on the other hand, impacts all 
energy terms and results in |AF^totai| ^ 7.403 MeV 
for setup (I). This difference is larger than for ^^^Pb. 
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L [fm] 

sim. 

resol. 

Ebind 

Etotal 

Ekin 

Eo 

El 

E 2 

E 3 E 4 Ec 

Sky3D 

24 

(I) 

1 fm 

-7.521 

-1790.095 

4493.687 

-19719.820 

315.821 

264.729 

11991.680 -90.798 954.640 

Sky3D 

24 

(I) 

0.5 fm 

-7.521 

-1790.056 

4493.707 

-19720.940 

315.855 

264.758 

11992.730 -90.799 954.640 

M-SHF 

200 

(I)* 

10 “'> 

-7.548 

-1796.439 

4509.010 

-19808.196 

316.196 

264.335 

12057.925 -91.182 955.474 

M-SHF 

200 

(I) 


-7.517 

-1789.036 

4493.844 

-19716.954 

315.797 

264.682 

11989.662 -90.682 954.613 

M-SHF 

200 

(I) 

10 “^ 

-7.517 

-1789.037 

4493.815 

-19716.783 

315.794 

264.672 

11989.534 -90.681 954.612 

Sky3D 

24 

(ii) 

1 fm 

-10.284 

-2447.609 

4524.698 

-19932.520 

322.707 

267.612 

12160.020 -92.395 302.305 

Sky3D 

24 

(ii) 

0.5 fm 

-10.284 

-2447.561 

4524.779 

-19933.900 

322.747 

267.675 

12161.240 -92.414 302.315 

M-SHF 

24 

(II)* 

10 "'^ 

-10.318 

-2455.747 

4538.262 

-19997.608 

322.283 

267.574 

12207.267 -94.137 300.612 

M-SHF 

24 

(II) 

10 "® 

-10.282 

-2448.297 

4523.168 

-19918.433 

322.321 

267.700 

12149.303 -92.444 300.090 

M-SHF 

24 

(II) 

10 "'^ 

-10.288 

-2448.444 

4523.097 

-19917.579 

322.290 

267.709 

12148.620 -92.475 299.893 


TABLE III: Parameters and energies for SkySD and M-SHF simulations of the nucleus. Table parameters are the same 
as in table 1^ The experimental binding energy is Ebind ~ 7.570 MeV [56] . 
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EIG. 4: Zoom of Eig|^b) showing the p, r and Ap prohles 
along X for the ^°®Pb simulation (I) with smoothing and e = 
10“^ (marked by a ★) and e — without smoothing. While 
the oscillation patterns of the three quantities are similar with 
and without smoothing, small differences are visible in the 
amplitudes. 
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FIG. 5: Evolution of the maximum error and binding energy 
with iterations of the ^°®Pb calculation. 


Figure [^b) shows a comparison between the x-profiles 
of p, r and Ap for the converged state with e = 10“^ 
and e = 10“^ whereas the latter contains smoothing. 
Small differences are present, especially in Ap. For 
better visualization we zoom in again and plot p, r and 
Ap for 0 fm < X < 8 fm, 0.145 fm < p < 0.18 fm and 
—0.025 fm < p < 0.025 fm for p, r and Ap, respectively 
in Fig. We can now see that the differences between 
simulations with and without smoothing are more 
pronounced than for ^^^Pb and we assume that this 
leads to the larger value of |AF^totai|- 

Since its initial shape is very different from the final 
configuration, has to undergo many shape changes 
during the iterations. On eight CPU nodes this results 
in simulations times on the order of days for about 4500 
iterations with a final truncation threshold of e = 10“^. 
The evolution of \og{S'ip) and F^bind are shown in Fig. 
1^ Unlike the ^^^Pb simulation, \og{52p) has now many 
local maxima and minima. These correspond to shape 
changes of the nucleus until the lowest energy state 
is found. The evolution of the binding energy follows 
a gradual decrease with some small fiuctuations. The 
latter seem to be a result of the interplay of Gaussian 
smoothing with the evolving wave functions. As soon 
as smoothing is removed around iteration 1500, the 
fiuctuations in the binding energy disappear and the 
latter jumps to higher values. The experimental value 
for ^^^U is F^exp = —7.570 MeV and thereby in good 
agreement with the simulations. 

As before, simulations of type (II) are initialized in 
the same way as simulations (I). For e = 10“^ and no 
Gaussian smoothing, both calculations setups evolve 
similarly up to iteration ^ 4400. At this point, as 
is shown in Fig(§ Sip of setup (II) stops to decrease 
and increases again. It reaches a maximum around 
iteration 8800 and then decreases slowly until iteration 
21000 where we stop the simulation. At this point, the 
configuration has not yet reached convergence according 
to our criteria, however, S'lp is very small, around 10“^ 
and is unlikely to grow again. Upon examining the 
cause of the increase in Sip we find a change in the 
orientation of the ^^^U nucleus. Figure IT shows its 
density iso-surface corresponding to p = 0.08 fm“^ at 
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distance [fm] distance [fm] 


FIG. 6: (a) Number density profiles of the along the x, y, and z-axis at convergence (subscript /) and the beginning of 
the iterations (subscript i). Subfigure (b) shows the number density, kinetic density r and laplacian of the number density Ap 
at iteration 2800. 



distance [fm] 


FIG. 7: Zoom of Fig|^b) showing the p, r and Ap profiles 
along X for the simulation (I) with smoothing and e = 

10“^ (marked by a 'k) and e = 10“^ without smoothing. 


iteration 4000 (left subfigure) and iteration 19400 (right 
subfigure). We can see that while the symmetry axis of 
the nucleus is initially parallel to the x-axis it changes 
to a diagonal of the simulation volume. This orientation 
seems to be energetically more favorable for periodic 
boundary conditions although the difference in binding 
energies between iteration 4000 and iteration 194000 
is very small, around |AF^bind| ^ 0.0106 MeV, as can 
be seen in Figj^ It is encouraging that our code can 
find the energetically more favorable state. However, 
the required timescales and iteration numbers are very 
large. Simulations with e = 10“^, with and without 
Gaussian smoothing, evolve similarly whereas they don’t 
reach convergence within 21000 iterations. Furthermore, 
their symmetry axes are not perfect diagonals of the 
simulation volume as is the case for the calculation with 




FIG. 8: Evolution of the maximum wavefunction change ^'0 
and binding energy Ebind with iterations for the nucleus. 



FIG. 9: Evolution of the maximum wavefunction change d'ljj 
and binding energy Ebind with iterations for the nucleus 
for simulation setup (II), e = 10~^ and no Gaussian smooth¬ 
ing. 
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FIG. 10: Density iso-surface of the nucleus for p — 

0.08 fm^ for simulation (II) with e = 10~^ and no Gaussian 
smoothing. The left figure shows the orientation of the nu¬ 
cleus in the simulation space at iteration 4000. The right 
figure shows the nucleus at iteration 19400. The black lines 
mark the symmetry axis. 


e = 10 The energies of all three calculations are 
given in table |m| Interestingly, different to simulation 
(I), the nucleus calculated with e = 10“^ in the 

MADNESS simulation with periodic be is more bound 
than for Sky3D. This is most likely due to the discussed 
rotation of the nucleus. Although the Sky3D simulation 
ran for about 30 000 iterations, the orientation of the 
symmetry axis stayed parallel to the x-axis. As 
before, the nuclear configurations in simulations (II) are 
more bound than in simulations (I) due to the presence 
of jellium. 

At this point, the findings regarding our M-SHF 
code can be summarized as follows: For small nuclei 
such as the code is fast and can determine the 

ground state within a few iterations to high precision, 
for both, small and large simulation volumes. Although 
both converged nuclei, and ^^^Pb, are spherically 
symmetric making harmonic oscillator states are a good 
initial guess, the ^^^Pb simulations requires significantly 
more computational time. In MADNESS, the latter 
scales with truncation threshold e which we adjust 
during the computation. Proper timing tests with a 
fixed truncation threshold will be done in the future 
and will give more detailed information on the scaling 
of the code. From our simulations, we see that 

the M-SHF code can find the ground state of nuclear 
configuration through several shape and orientation 
changes and is in agreement with experimental binding 
energies and other SHF simulations. With that, we turn 
to the study of nuclear pasta. 


VI. NUCLEAR PASTA FROM MOLECULAR 
DYNAMICS 

A. Without spin-orbit contributions 

Numerical Skyrme Hartree-Fock studies of nuclear 
pasta phases are usually performed by initializing single 
particle states as plane waves or random positioned 
Gaussians. Ideally, the nucleon wave functions then 
converge into the ground state configuration. However, 
as matter frustration allows for many different local 
energy minima and thereby pasta shapes, matter can 
easily become trapped in a quasi-ground state. To 
facilitate the search for the true ground state at a given 
density, proton fraction and temperature, restrictions 
can be placed on the symmetry and shape of the nuclear 
configuration m- This leads to a faster convergence of 
the latter and, by changing the symmetry assumptions, 
allows to scan through different pasta shapes. The 
ground state can then be identified as the one with the 
lowest total energy. The drawback of this approach is 
that the final shapes are somewhat predetermined by 
the specific assumptions and it might be difficult to 
explore new geometries. 

Molecular dynamics simulations start out with randomly 
placed nucleons that are evolved over many iterations. 
The only restrictions for such methods are the imposed 
boundary conditions and simulation space dimensions. 
Different variations of the latter can be tested to ensure 
that the ground state is independent of the simulation 
space setup. Especially large molecular dynamics 
calculations with > 10^ nucleons and box lengths of 
L > 100 fm minimize the effects of the simulation space 
geometry. The resulting pasta configurations can be tra¬ 
ditional or novel shapes as found in [32]. However, MD 
simulations often do not contain quantum mechanical 
features and their results should be cross-checked with 
self-consistent calculations such as Hartree-Fock. In 
this work, we use the converged pasta configuration of 
a simulation with the the Indiana University Molecular 
Dynamics code lUMD [131 El 1221 1221 l58] and explore 
its evolution as we iterate the single particle states with 
the M-SHF code. 

Our starting point is the so-called waffie phase [32] - an 
intermediate state between the lasagna (plate) and the 
spaghetti (rod) configurations. It consists of plates with 
a lattice of periodic holes whereas two neighboring plates 
are displaced by half of the lattice spacing. The lUMD 
calculation was performed using periodic boundary 
conditions with 490 neutron and 210 proton particles in 
a simulation box with length L = 24 fm. The average 
density is p = 0.05 fm~^ with a proton fraction Yp = 0.3. 
The converged nucleon positions are shown in Fig[TT] 
with blue (dark) spheres symbolizing neutrons and grey 
(light) spheres protons. Note that the temperature 
of the MD simulation was 1 MeV while our M-SHF 
calculations are performed at zero temperature. Typical 
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FIG. 11: Positions of neutrons (blue) and protons (grey) in the converged lUMD pasta simulation. Subfigures show different 
orientations of the simulations space. Sphere sizes are for visualization only. 



FIG. 12: Iso-surfaces of initial total baryon number density p(f) of M-SHF. Nucleon wavefunctions are Gaussians folded around 
coordinates from lUMD shown in Fig[^ Subfigures correspond to orientations of the simulations space as in Fig[^ See text 
for details. 


pasta calculations are performed either at supernova 
conditions - finite temperature and ^ 0.2 — 0.4 - or 
neutron star conditions with T ^ 0 and Yp < 0.1. The 
configuration that we are describing, i.e. T = 0 MeV 
and Yp = 0.3, does therefore not exist in neutron stars 
or supernova environments. However, the aim of this 
paper is to explore the stability of the lUMD pasta 
configuration once quantum mechanical features are 
added and to provide a proof-of-principle study for 
future finite temperature M-SHF simulations that will 
apply neutron star or supernovae conditions. 

For the initialization, each single particle state is given 
by a sum of 27 3D Gaussians with a = 3.0 fm. The 
Gaussian are centered at the nucleon coordinates from 
lUMD and their closest images due to the reflective 
boundary conditions. Figure shows the resulting 
density iso-surfaces for 0.08 fm“^^< p < 0.16 fm“^ and 
different orientations of the simulation space. We can 
identify two plates, each with one hole. The latter 
are displaced relative to each other so that a hole is 
aligned with a denser region in the neighboring plate. 
The M-SHF iterations are performed as previously 
discussed for setup (H). Here, we do not include either 
the spin-orbit potential Uso or the spin-orbit density J 
and will test their effects in the next section. 

As for the nuclear ground states in the previous dis¬ 
cussion, we plot the evolution of the maximum wave 



Iterations 

FIG. 13: Evolution of the maximum wavefunction change S'ljj 
and binding energy Ebind with iterations for the MD M-SHF 
pasta simulation 


function change and binding energy F^bind Fig.[T^ 
with vertical lines indicating the reduction of e. The 
simulation was stopped after ^ 8000 iterations when 
^ 10“^. It ran for about two weeks on 24-30 nodes. 

From the maxima and minima in we see that 
the pasta shapes underwent some small changes in the 
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FIG. 14: Iso-surfaces of the total baryon density p(r) of the converged M-SHF simulation at iteration 8000. Subhgures are as 
in Figp^. 



L [fm] 

sim. 

resol. 

Ebind 

Etotal 

Ekin 

Eq 

El 

E 2 

E 3 

E4 Ec 

M-SHF 

24 

(H) 

10 “' 

-9.059 

-6341.536 

12997.724 

-51349.350 

529.618 

449.606 

30803.440 

- 227.427 

Sky3D 

24 

(11) 

1 fm 

-9.083 

-6358.075 

13019.705 

-51434.630 

525.969 

451.538 

30853.410 

- 225.933 


TABLE IV: Parameters and energies for the nuclear pasta simulation with M-SHF and SkySD. The table setup is as in table [I] 


beginning of the simulation for iterations < 4600. Since 
the reduction of the truncation threshold from e = 10“^ 
to e = 10“^ occurs very early in the simulation, we 
decided to apply Gaussian smoothing until the second 
reduction of e to 10“^. This explains the jump in 
binding energy around iteration 2800. The binding 
energy then converges to a value of F^bind ^ —9.059 MeV. 


The different energy components are given in table IV 


together with the results of the Sky3D pasta simulation. 
The latter also initialized the wavefunctions also via 
Gaussians folded around the MD nucleon coordinates. 
The Sky3D calculation converged after ca. 32000 
iterations. The difference in total binding energies for 
Sky3D and M-SHF is |AF^totai| ^ 16.539 MeV which is 
about ^ 2.60 x 10“^|F^totai|- 

Figure shows the density iso-surfaces, again for 
0.08 fm“^ < p < 0.16 fm“^ in the final M-SHF pasta 
configuration at iteration 8000. As before, subfigures 
correspond to different box orientations. The general 
shape of the waffle phase is very similar to Fig.p!^ 
Main differences are the broadening of regions with 
p > 0.08fm~^, while the holes in both plates are smaller. 
The variations of correspond to these changes, 
whereas the smooth decrease of for iterations > 4600 
implies that the ground state configuration at iteration 
8000 is already reached at this point. Figure [T^a) 
compares the number density profiles of the pasta 
configuration at iteration 4600 and 8000. The differences 
are very small which also applies to the kinetic densities 
and the laplacian of the densities along the x-axis, as 
shown in Fig|T^b). For future studies we should there¬ 
fore consider to either modify the convergence criteria, 
consider a configuration to be stable at an earlier point 
in the simulation, or change the iteration procedure so 
that the convergence criteria are met faster once the 
pasta shape does not undergo significant changes. 



-10 -5 0 5 10 


distance [fm] 



FIG. 15: (a) Number density prohles of the nuclear pasta 
conhguration along the x, y, and z axis at iteration 4600 and 
in the hnal, i.e. converged, state at iteration 8000. (b) Gom- 
parison of the x-prohles of the total number density p, kinetic 
density r and laplacian of the number density Ap in the h- 
nal, i.e. converged, state of the nuclear pasta and at iteration 
4600. 
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FIG. 16: Evolution of the maximum wavefunction change d'ljj 
and binding energy Ebind with iterations for the MD M-SHF 
pasta simulation 


B. With spin-orbit contributions 

Although, the spin-orbit contribution is not expected 
to be a large part of the total nuclear binding energy, 
it is important for the reproduction of energy levels and 
magic numbers. It has also been found that the spin- 
orbit potential can impact the deformation of nuclei [59] , 
however, the effect on nuclear pasta phases has not been 
fully explored yet. 

In this section we study the impact of the spin-orbit po¬ 
tential on the waffle pasta phase. We perform the calcu¬ 
lations with Sky3D and M-SHF. For the latter, the initial 
configuration is taken from the previous section at iter¬ 
ation ^ 2000 when d'ljj ^ 10“^ (see Figj^. For Sky3D, 
the starting point is the converged MD state. Here, as 
in the previous study, Gaussians are folded around the 
nucleon coordinates and evolved. However, now the nu¬ 
clear potential contains the spin-orbit terms from the 
very beginning. For M-SHF, FigfT6| shows the evolution 
of 5^1; and F^bind as functions of iteration number. Verti¬ 
cal dashed lines indicate the reduction of the truncation 
threshold. Note that the reduction of e from 10“^ to 
10“^ happens quite late in the simulation. Due to speed 
we evolved the configuration to lower with 10“^ and 
reduced the truncation threshold only at the very end of 
the simulation. However, the switch from e = 10“^ to 
10“^ and consequent evolution over ca. 400 iterations 
did not cause any changes in the pasta configuration or 
its convergence. 

A large jump can be seen in as well as a step¬ 
like increase in |F^bind| as soon as we add the spin-orbit 
terms. The second jump around iteration 3500 when we 
reduce the truncation threshold from e = 10“^ to 10“^ 
is due to the removal of Gaussian blurring. The simu¬ 
lation converges until iteration ^ 7500 where starts 


to increase again. The different minima and maxima be¬ 
tween iterations 10000 - 18000 indicate possible shape 
changes. Eventually, starts to continuously decrease 
and reaches convergence with r\j 10 ^ for e = 10 
In total, the simulation requires about 21000 iteration 
steps whereas the Sky3D calculation converges already 
after iteration 6000. It is not clear whether adding the 
spin-orbit potential from the very beginning would also 
lead to a quicker convergence for M-SHF. This has to be 
explored in the future. 

Table |V| compares the final energies of the two simula¬ 
tions. The energy contributions are similar. As before, 
we find that the absolute of the binding energy is smaller 
for M-SHF than for Sky3D. The difference is about 
lA^totail ^ 5 MeV which is only ^ 7.77 x 10“"^|F;totai|- 
With spin-orbit, the waffle phase is more bound, by about 
- 82.72 MeV and - 94.26 MeV for Sky3D and M-SHF, 
respectively, which is ^ 1.3 — 1.5% of the total energy. 


Figures [T7| and [18| show the converged pasta for M-SHF 
without and with spin-orbit terms and the Sky3D sim¬ 
ulation with spin-orbit. Although the general shape is 
the same, we can find subtle differences. Without spin- 
orbit, the size of the holes in the top and bottom plates 
seem very similar. When the spin-orbit contributions 
are added, one hole shrinks while the other one becomes 
larger. This evolution corresponds to the shape changes 
that were indicated in Fig.[T6| between iteration 10000 - 
18000. Interestingly, the M-SHF and Sky3D simulations 
both show this effect despite the different initializations. 
For M-SHF, the small hole is in the top plate while the 
large one is in the bottom plate. For Sky3D the situation 
is reversed. However, due to the periodic boundary con¬ 
ditions the order is not important and the pasta phase 
should consist of a lattice of alternating small and large 
holes. The question is of course, whether the same struc¬ 
ture would be found in simulations of a larger volumes 
or with different spin-orbit potentials. More systematic 
studies have to be performed in the future. At present, 
we conclude that the inclusion of the spin-orbit contribu¬ 
tion in the Sv-bas nuclear potential modifies features of 
the waffle phase but does not lead to its disappearance. 


VII. SUMMARY 

In this work, we introduce and discuss calculations 
of nuclear matter via Skyrme Hartree-Fock calculations 
with the Multi-resolution ADaptive Numerical Environ¬ 
ment for Scientific Simulations (MADNESS). To verify 
and benchmark the code, we perform calculations of nu¬ 
clear ground states and find good agreement with the 
established Skyrme Hartree-Fock code Sky3D and ex¬ 
perimental binding energies. While calculations for light 
nuclei seem to be very fast, the scaling of the code with 
number of nucleons needs improvement for future stud¬ 
ies. We test our code for large boxes and free bound- 
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L [fm] 

sim. 

resol. 

Ebind 

Etotal 

Ekin Eq 

El 

E 2 

E3 E4 Ec 

M-SHF 

24 

(n) 

10 “' 

-9.194 

-6435.792 

13327.675 -52750.085 

550.097 

486.129 

31884.606 -168.597 234.383 

Sky3D 

24 

(ii) 

1 fm 

-9.201 

-6440.795 

13329.782 -52767.190 

549.622 

487.920 

31897.320 -171.994 233.754 


TABLE V: Parameters and energies for the nuclear pasta simulation with M-SHF. The table setup is as in table 



10 5 0 -5 -10 

X rfml 



10 5 0 -5 -10 

X [fml 



X rfml 


FIG. 17: Iso-surfaces of p{r) as in Fig[T4] The orientation corresponds to the top plate of the converged waffle phase. Subfigure 
(a) shows the result of the M-SHF simulation without spin-orbit and (b) with spin-orbit interactions. Subfigure (c) is the 
converged SkySD simulation with spin-orbit. 


ary conditions and small boxes with periodic boundary 
conditions. For nuclear pasta simulations, we explore 
a configuration in a 24 fm box with periodic boundary 
conditions and 700 nucleons with a proton fraction of 0.3 
and an average density of p = 0.05 fm“^. The initial¬ 
ization of the simulation is done using the output of a 
converged simulation by the Indiana University Molecu¬ 
lar Dynamics code. The corresponding shape is the waf¬ 
fle phase. For a calculation without spin-orbit terms, we 
find that the simulation fulfills our convergence criterium 
after 8000 iterations whereas the configuration and bind¬ 
ing energies do not change significantly after iteration 
^ 4000. Furthermore, the final shape of the nuclear 
configuration does not differ significantly from the initial 
MD state indicating that the waffle phase is stable even 
when quantum mechanical effects are considered. When 
adding spin-orbit nuclear potential terms to a partially 
converged calculation, we find small shape changes which 
push the convergence to iteration ^ 20000. The shape of 
the waffle phase has small but visible differences in com¬ 
parison to the calculation without spin-orbit. However, 
the phase remains in the waffle geometry. Similar results 
are also found with the Sky3D calculation. 
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